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Abstract 

The  time  independent  electric  potential  due  to  water  and  a  lithium  ion  near  a  charged  metal 
surface  is  calculated  by  space  and  ensemble  averaging  of  trajectories  generated  by  a  molecular 
dynamics  simulation.  Since  the  cation  does  not  contact  adsorb  variations  in  the  electric 
potential  near  the  metal  surface  are  due  to  water  oriented  in  the  electric  field  of  the  charged 
surface.  The  potential  is  decomposed  into  separate  contributions  from  monopoles  (from  the 
ions),  and  dipoles,  quadrupoles  and  octopoles  (from  the  water  molecules).  At  distances 
greater  than  about  0.5  nm  from  the  electrode  (2  -  3  water  molecules)  the  potential  is  'flat' 
with  the  quadrupole  contributing  most  due  to  a  near  cancellation  of  the  ion  and  water  dipole 
components.  Approaching  the  surface  weak  features  are  encountered  due  to  water  packing 
and  then  a  big  oscillation  due  to  water  oriented  in  a  layer  next  to  the  electrode.  None  of 
these  effects  are  described  in  theories  that  approximate  water  as  a  continuum  fluid. 
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Abstract 

The  time  independent  electric  potential  due  to  water  and  a  lithium  ion  near  a  charged  metal 
surface  is  calculated  by  space  and  ensemble  averaging  of  trajectories  generated  by  a  molecular 
dynamics  simulation.  Since  the  cation  does  not  contact  adsorb  variations  in  the  electric  potential 
near  the  metal  surface  are  due  to  water  oriented  in  the  electric  field  of  the  charged  surface.  The 
potential  is  decomposed  into  separate  contributions  from  monopoles  (from  the  ions),  and  dipoles, 
quadrupoles  and  octopoles  (from  the  water  molecules).  At  distances  greater  than  about  0.5  nm 
from  the  electrode  (2  -  3  water  molecules)  the  potential  is  'flat'  with  the  quadrupole  contributing 
most  due  to  a  near  cancellation  of  the  ion  and  water  dipole  components.  Approaching  the  sur¬ 
face  weak  features  are  encountered  due  to  water  packing  and  then  a  big  oscillation  due  to  water 
oriented  in  a  layer  next  to  the  electrode.  None  of  these  effects  are  described  in  theories  that 
approximate  water  as  a  continuum  fluid. 
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I.  INTRODUCTION 

In  this  paper  we  describe  the  calculation  of  the  electric  potenUal  across  an  electrified  interface 
using  a  molecular  dynamics  simulation  performed  at  constant  N,V,T.  The  simulation  cell  con¬ 
sists  of  one  lithium  ion  and  157  water  molecules.  Lithium  was  chosen  because  it  does  not  lose 
solvation  and  contact  adsorb  on  the  surface  in  the  range  of  charge  densities  used  (-0.288e/nni^) 
in  the  calculation.  This  makes  the  interpretation  of  structure  near  the  electrode  simpler.  We 
resolve  the  potential  into  explicit  contributions  from  the  ions  and  the  multipole  moments  of  the 
water  molecules.  The  variation  in  the  potentials  with  distance  from  the  electrode  shows  structure 
which  is  related  the  distribution  of  the  ion  and  the  distribution  and  orientation  of  the  water 
molecules.  The  features  in  the  electric  potential  due  to  water  are  new  and  not  found  in  theories 
that  treat  the  solvent  as  a  dielectric  continuum.  The  explicit  calculation  of  quadrupole  and 
octopole  potentials  provides  insight  not  available  from  theories  that  model  the  solvent  as  a  point 
dipole  in  a  sphere. 

Figure  1  upper  part  shows  schematically  the  electric  potential  of  the  system  when  water  is  as¬ 
sumed  to  be  a  simple  dielectric  continuum.  The  lower  part  of  Figure  1  is  a  cartoon  that  provides 
on  the  molecular  scale  a  picture  of  the  double  layer  near  a  flat  charged  metal  surface  when  there 
are  no  contact  adsorbing  ions.  Similar  pictures  can  be  found  in  many  textbooks  and  review  ar¬ 
ticles.  The  water  next  to  die  electrode  is  shown  as  an  oriented  monolayer.  If  this  first  layer  of 
water  is  not  displaced  by  the  cation  at  its  distance  of  closest  approach  then  the  hydrated  cation 
is  two  water  molecules  away  from  the  surface.  This  is  a  characteristic  feature  of  the  model  of 
Bockris,  Devanathan  and  Miillerl-^.  In  Grahame's  model'^  the  ions  primary  solvation  shell  can 
contact  the  electrode.  The  plane  of  closest  approach  of  the  positive  ion  is  the  outer  Helmholtz 
plane  (OHP).  The  diffuse  region  in  this  traditional  picture  starts  two  solvent  molecules  from  a 
flat  electrode  surface  and  stretches  out  many  nanometers  into  the  bulk  electrolyte.  The  Gouy 
and  Chapman  modeP  provided  the  first  quantitative  description  of  the  diffuse  layer. 

This  paper  consists  of  three  sections.  Section  II  contains  a  brief  description  of  model  and 
methods.  Two  complementary  methods  were  used  to  calculate  the  electric  potential.  In  the  last 
section  III,  we  first  briefly  describe  the  distribution  profiles  for  ions  and  water  molecules  ob¬ 
tained  from  the  molecular  dynamics  simulation,  and  then  the  results  of  the  electric  potential 
calculations.  These  include  the  results  of  decomposing  the  potential  into  multipole  components 
and  a  summary  of  the  effect  of  varying  the  width  of  tlie  spatial  box  used  in  making  the  space 
averages.  It  is  shown  how  the  water  quadrupole  and  octopole  multipole  structure  in  the  potential 
is  washed  out  when  the  box  dimension  parameter  approximates  the  size  of  a  water  molecule. 

II.  MODEL  AND  METHODS. 

A.  The  screened  electrode 

Since  the  electrolyte  solution  is  a  conductor  all  charged  surfaces  immersed  in  the  electrolyte  are 
screened  by  displacement  of  ions  towards  the  surface  in  such  a  way  as  to  shield  the  bulk  solution 
from  the  field  of  the  surface  charge.  The  simplest  treatment  of  this  phenomena  is  the  Gouy- 
Chapman  theory^  which  treats  ions  as  point  charges  and  all  the  solvent  molecules  as  if  they 
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comprised  a  dielectric  fluid.  According  to  this  theory  the  concentration  of  monovalent  ions  has 
fallen  to  e'*  of  its  'surface'  value  at  a  distance  d  from  the  electrode  given  by 


pii 

y  \  Sne^ni,  / 

The  inverse  d'^  is  the  Debye-Hiickel  screening  constant.  Here  e  is  the  macroscopic  dielectric 
constant  of  the  solvent,  and  nj^  is  the  concentration  of  the  ions  in  the  bulk.  Typical  values  of  d 
are:  3.1  nm  for  0.01  M,  and  0.96  run  for  0.1  M.  The  derivation  of  this  formula  does  not  hold 
high  salt  concentrations,  though  it  can  provide  some  guidance.  For  example  at  0.3  M  we  have 
d  =  0.55  nm  which  is  only  about  twice  the  size  of  a  water  molecule  or  ion.  We  conclude  that 
most  of  the  surface  charge  screening  takes  place  close  to  the  electrode  and  that  a  simulation  cell 
with  size  L  =  1.862  nm  should  be  just  large  enough  to  capture  most  of  the  physics. 


B.  Immersed  electrode  model. 

The  model  consists  of  a  layer  of  electrolyte  between  two  walls.  The  gap  between  the  walls  is 
the  simulation  cell  edge  length  L  =  1.862  nm.  The  cell  contains  157  model  ST2  water  molecules 
and  one  lithium  ion.  The  cell  is  periodically  replicated  in  the  two  dimensions  xy  parallel  to  the 
metal  surface,  but  not  perpendicular  to  the  surface.  The  wall  on  the  left  (z  =  -L/2)  carries  no 
charge,  it  is  a  simple  restraining  wall  to  hold  the  electrolyte  in  place.  The  wall  on  the  right  (z 
=  L/2)  is  the  metal  with  the  plane  z  =  L/2  the  electrostatic  image  plane.  The  charge  on  the  metal 
equals  the  image  charge  of  the  lithium  ion.  The  total  charge  of  the  system  is  zero  because  the 
water  molecules  are  individually  neutral  and  the  ion  and  its  electrostatic  image  are  opposite  in 
charge.  The  essential  feature  of  the  immersed  electrode  model  is  that  the  charge  in  solution 
exactly  equals  the  charge  on  the  metal,  and  the  charge  on  the  metal  is  the  net  electroststic  image 
charge  (all  the  water  molecules  have  electrostatic  images  but  their  net  charge  is  zero). 

C.  Model  for  water  and  ions. 

In  all  the  calculations  reported  here  we  use  the  parameters  of  the  Stillinger  ®  ST2  water  model 
and  the  interaction  parameters  for  alkali  metal  ions  and  water  developed  by  Heinzinger  and  co¬ 
workers  The  ST2  water  molecule  model  consists  of  a  central  oxygen  atom  (0_ST2  or  O  for 
short)  surrounded  by  two  hydrogen  atoms  (H_st2  or  H  for  short)  and  two  massless  point  charges 
(PC_ST2  or  PC  for  short)  in  a  rigid  tetrahedral  arrangement.  The  0-H  and  0-PC  bond  lengths 
were  0.10  nm  and  0.08  nm  respectively.  This  small  difference  in  bond  lengths  means  that  the 
water_st2  model  and  its  electrostatic  image  behave  similarly.  The  only  Lennard- Jones  'atom'  in 
the  ST2  model  is  the  oxygen  atom.  The  hydrogen  H_ST2  and  point  charges  PC_ST2  interact 
with  their  surroundings  (i.e.  other  atoms  and  surfaces)  only  via  Coulomb  interactions.  Their 
charges  are  qf^=0.23570lel  and  qpc=-qH-  "Th^  O  atom  has  zero  charge.  The  lithium  ion  was 
treated  as  a  non-polarizable  Lennard-Jones  atom  with  point  mass  and  charge.  For  example  the 
(e,a)  pairs  are  (0.3164,  0.3100)  and  (0.1490,  0.2370)  for  0_ST2  and  Li  ion  respectively.  The 


units  are  e  in  kJAnole  and  o  in  run  .  The  Lorentz-Berthelot  combining  rules  were  enforced  for 
unlike  species,  namely;  Eab  =  {eMeBs)^  and  Oab  =  ¥2(0^  +  Gbb).  The  Stillinger  switching  function 
was  used  to  modify  the  coulomb  interactions  between  water  molecules  at  close  range,  ends  All 
molecule-molecule  Lennard- Jones  type  interactions  were  cut-off  in  a  smooth  fashion  at  a  mo¬ 
lecular  separation  R  =  0.68  nm. 

D.  Interaction  between  water  and  ions  and  the  walls. 

The  metal  was  represented  by  two  linearly  superimposed  potentials.  Pauli  repulsion  and 
dispersive  attractive  interactions  were  represented  by  a  9-3  potential,  and  the  interaction  with  the 
conduction  electrons  by  an  electroststic  image  potential.  In  the  calculations  described  here  the 
image  plane  and  origin  plane  of  the  9-3  potential  were  chosen  coincident.  This  was  tantamount 
to  choosing  the  image  plane  and  the  nuclear  plane  of  the  metal  surface  to  be  same  plane.  This 
is  acceptable  in  our  scheme  because  the  Lennard-Jones  core  parameters  a  are  all  large  and  the 
'thickness'  of  the  repulsive  wall  is  also  large  (ca.  0.247  nm).  The  atom-surface  interaction  pa¬ 
rameters  describing  interaction  with  nonconduction  electrons  were  chosen  to  be  the  same  as 
those  used  by  Lee  at  al  *,  The  well  depth  is  approximately  equal  to  thermal  energy  'kT'  or  about 
2.4  kJ/mole. 


G.  Electrostatics 

During  the  molecular  dynamics  simulations  the  fast  multipole  method  developed  by  Greengard 
and  Rokhlin  was  used  to  evaluate  all  long  range  sums  of  electrostatic  interactions  without 
truncation.  The  charge  on  the  electrode  is  the  image  charge  -lei.  At  any  instant  during  the 
molecular  dynamics  run,  this  charge  is  not  uniformly  distributed  across  the  electrode  but  local¬ 
ized  on  the  surface  in  such  a  way  as  to  produce  the  same  electric  field  and  potential  as  the 
electrostatic  image  of  the  lithium  ion  and  all  the  water  molecules.  The  field  acting  on  the  lithium 
ion  comes  from  all  the  water  molecules  in  the  cell,  all  water  and  ions  in  the  xy  periodical 
replicant  cells,  and  all  of  the  electrostatic  images  of  the  contents  of  all  cells  in  the  image  plane 
of  the  metal.  We  emphasize  again  that  in  this  calculation  as  in  al!  the  others  described  in  this 
paper  no  electrostatic  interaction  is  truncated. 

There  are  two  useful  ways  to  calculate  the  average  time  independent  electric  potentials  from 
stored  molecular  dynamics  trajectories.  The  reason  for  doing  this  is  to  make  comparisons  with 
averaged  potentials  calculated  by  other  theories,  e.g.,  Gouy-Chapman.  The  first  method  (the 
atom  approach)  uses  the  charge  or  partial  charge  on  each  atom  and  the  three  space  coordinates 
of  each  atom  recorded  at  regular  lime  separations  (usually  1  ps)  as  input  into  the  generation  of 
a  charged  source  distribution.  The  source  term  for  Maxwell  s  equations  is  the  charge  density 
function  in  vacuum.  Without  any  local  spatial  averaging  this  point  of  view  is  similar  to  that 
followed  by  Wilson,  Pohorille,  and  Pratti^’  The  second  way  (the  molecule  method)  views 
the  dynamical  system  as  a  collection  of  molecules  with  internal  electrical  structure.  These  are 
inherently  more  complex  objects  than  atoms.  To  specify  the  electrical  properties  we  need  the 
space  coordinates,  orientations  and  electrostatic  multipole  moments  of  each  ion  and  molecule  in 
the  system.  The  approach  we  follow  was  described  by  Russakoffl^.  This  derivation  of  a  set 
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of  macroscopic  Maxwell  equations  from  the  equations  of  charged  particles  in  vacuum  is  sum¬ 
marized  in  a  number  of  standard  texts^^.  We  note  that  though  the  dynamics  are  uniquely  defined 
by  the  models  we  use.  the  time  independent  average  fields  and  potentials  that  are  calculated  are 
not  unique  because  we  can  broaden  the  point  charge  distributions  inside  the  ions  and  molecules 
in  a  spherically  symmetric  way  and  so  long  as  the  broaden  distributions  stay  well  inside  the 
Lennard-Jones  spheres  (so  that  distributions  on  different  molecules  do  not  overlap)  the  the  dy¬ 
namics  are  unchanged.  This  means  we  can  smooth  some  of  the  averages  within  limits  set  by 
geometry  of  the  models  and  the  energetics  of  collisions.  This  point  has  already  been  pointed 
out  by  Wilson,  Pohorille  and  Pratt^^’  in  their  interesting  discussion  of  the  properties  of  the 
vacuum-water  surface. 

The  atom  method  uses  the  distribution  of  point  charges  in  vacuum.  Consider  the  set  of  point 
charges  without  regard  for  whether  they  originate  from  neutral  water  molecules  or  charged  ions. 
In  this  case  the  source  term  for  Maxwell  fields  in  vacuum  is  the  microscopic  charge  density 


where 

=  p„^,„Xx,y,05(z  -  VlL), 


Paion:s(r,t)=  r/t)], 

y=i 


[2.2] 

[2.3] 


[2.4] 


Here  Kicms  is  the  number  of  atoms  in  the  simulation  cell  and  r;  is  the  position  of  the  j-th  atom. 
The  surface  charge  density  on  the  metal  p™,.;(x,y,0  is  time  dependent  because  it  depends  on  the 
position  of  the  atomic  charges.  All  the  charge  in  the  system  is  described  by  p;.  We  can  subject 
this  charge  density  to  local  space  averaging  using  a  test  function.  Let/(r)  be  a  real  positive 
function  localized  around  r  =  0.  We  define  the  local  spatial  average  of  F(r,  t)  by 


<F(r,  t)>  =  jdr'fir')F{r-r',  t). 


[2.5] 


The  time  average  can  be  replaced  by  an  average  over  configurations  at  different  times 


F(r)  = 


lim 

r  — >  oo 


^conftgs 


N. 


configs 


X  t,)>. 


i=  1 


[2.6] 
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The  system  has  translational  invariance  in  the  xy  plane  therefore  F{r)  is  only  a  function  of  z. 
Consequently  we  consider  only  test  functions  f{z)  like  the  localized  one  dimensional  Gaussian 
function 


Az)  =  (KgY^e-^<  f2.7] 

We  will  subsequently  refer  to  g  as  the  gaussian  width  or  'bin'  width.  The  metal  surface  charge 
density  is  replaced  by  the  averaged  image  charge  density  which  is  a  constant.  The 

metal  image  charge  is  denoted  by  After  the  configurational  averages  are  performed  we 

get  the  z  dependent  charge  density 


P/(z)  =  p„eM/6(z  -  V2L)  +  patoms(^) 


where  after  substituting  explicitly  for  the  test  function  we  get 


1 


N, 


configs 


I 

/  =  1 


^gy/[z-Zy(t,)]. 

;=i 


[2.9] 


The  electric  field  and  potential  are  given  by  integrating  the  vacuum  Maxwell's  equation  over  the 
remaining  space  variable  z 

£y/z)  =  -^j  dz'piiz')  ,  dz%Az').  [2.10] 

The  second  approach  treats  the  system  as  a  collection  of  molecules  (labels  n,  positions  r„) 
composed  of  atoms  (label  b,  position  charge  q„t).  The  charge  density  of  the  system  is  still 
the  same  as  p/(r,  /)  but  with  the  atoms  collected  into  molecular  groups.  For  the  molecule  method 
we  write  the  charge  density  as  follows: 


p,/(r,0  =  p,„<,to/(-r,y,r)5(z  -  ViL)  +  p,„o/(r,0 


[2.11] 


where  p,«a/  is  the  same  charge  density  on  the  metal  surface  as  in  the  atomic  method,  and  the 
charge  density  from  the  molecules  is 


=  ^p„(r,0 

H  =  1 


^X'^«i,5[r-r„(/)-r„fc(?)], 

n=  \b~\ 


[2.12] 
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Here  N„^  is  the  total  number  of  molecules  in  the  system.  The  number  of  atoms  in  the  n  th 
molecule  is  AL,  and  r„6  is  the  position  of  the  charge  <?„(,  measured  from  the  center  r„  of  the  n  th 
molecule.  Note  that  r„  and  r„b  are  time  dependent. 

Next  we  perform  explicit  local  spatial  averaging  with  a  test  function  and  take  the  ensemble  av¬ 
erage.  These  steps  are  formally  the  same  as  described  for  the  atom  method.  Then  for  each 
molecule  (label  n)  we  take  the  average  with  respect  to  the  test  function  y(z),  and  make  a  Taylor 
expansion  of  atomic  coordinates  r„i,  relative  to  the  molecular  centers  r„.  This  yields  the  fol¬ 
lowing  expression  for  the  averaged  charge  density  by  the  molecule  method 


2  3 

=  P(^)  -  ^  QM  -  -f  •  •  •  [2. 1 3] 

dz  dz 

where 


p(z)  = 


1 


^configs  N^l 


,-=l  n,  =  l 


6=1 


[2.14] 


1 


N, 


configs 


S 

!=  1 


m~l  6=1 


[2.15] 


^configs  Nmol  ^mb 

Qzziz)=-^ -  y  VlY  A2-2miti)]Tq,„bZmbitiKbifi) 

‘'‘configs  ^ 


[2.16] 


Nconfigs  N^iy 

Ozzzi^)  =  -1^ -  y  4“  y  ]y^qmb^mbU,Kb(q)^,nbiti)  [2.17] 

Nco^.gs  ^  6 

Here  the  charge  density  in  tlie  second  method  has  been  resolved  into  contributions  from 
monopoles  (ions  only  in  this  paper),  dipoles,  quadrupoles,  octopoles,  and  higher  order  terms  (all 
from  the  solvent  in  this  paper). 
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III.  SINGLE  METAL  ION 

In  this  section  we  first  describe  briefly  the  probability  density  distribution  calculated  from  mo¬ 
lecular  dynamics  trajectories.  Then  we  discuss  the  contibutions  that  the  electrostatic  multipoles 
make  to  the  total  electric  potential.  Finally  we  discuss  the  sensitivity  of  the  potential  and  its 
components  to  the  width  of  the  test  function  used  in  the  calculation. 

The  simulation  was  mn  for  2000  ps,  the  first  100  ps  of  which  were  used  to  equilibrate  the  sys¬ 
tem.  Figure  2  shows  the  probability  density  profiles  for  the  atoms,  ions  and  water  averaged  over 
the  xy  plane.  Also  shovwi  at  top  in  Figure  2  is  the  electrical  potential  calculated  by  the  atom 
method  using  a  gaussian  test  function  with  a  width  of  g  =  0.03  nm.  It  is  shown  with  the  prob¬ 
ability  distributions  to  highlight  features  due  to  water  the  majority  species.  Note  that  in  Figure 
2  the  scale  for  water  and  its  components  differ  from  the  lithium  ion  by  a  factor  of  fifty. 

In  Figure  2  the  Li"^  ion  mass  center  maps  out  a  diffuse-like  region  between  -0.6  and  0.4  nm. 
The  small  ionic  radius  of  Li"^  ensures  that  its  primary  solvation  shell  is  tightly  bound  making  it 
the  core  of  a  much  larger  composite  object.  Only  in  a  particularly  high  surface  electric  field 
can  this  object  dissociate  and  permit  the  LF  cation  to  contact  adsorb  on  the  electrode.  The  ap¬ 
parent  asymmetry  of  the  distribution  may  be  due  to  the  limited  width  of  the  cell  and  conse¬ 
quently  no  significance  is  attached  to  its  shape  other  than  it  is  diffuse  in  nature.  The  ion  can 
be  pulled  closer  to  the  metal  by  applying  an  external  uncompensated  electric  field.  Uncom¬ 
pensated  field  cannot  be  screened  completely  by  ions  present  in  the  solution.  This  is  convenient 
way  to  check  the  sensitivity  of  an  ion  to  adsorption  in  higher  fields,  without  having  to  resort  to 
a  separate  lengthy  computer  simulation.  When  the  surface  field  was  roughly  equal  to  that 
produced  by  doubling  the  negative  charge  on  the  metal,  the  lithium  ion  did  not  contact 
adsorb*^’  whereas  all  the  larger  cat  ions  (Na''^,  Rb''^,  and  Cs"'^)  did  contact  adsorb. 

To  check  for  large  effects  due  to  correlated  motion  between  two  Li"^  ions  in  the  same  simulation 
cell  we  performed  the  calculations  for  several  hundred  picoseconds  with  a  system  about  four 
times  larger  (two  lithium  ions  and  598  water  molecules).  The  water  density  profiles  when  ap¬ 
propriately  scaled  appeared  almost  exactly  the  same,  and  the  lithium  probability  p(z)  profile  was 
again  spread  across  the  cell  and  formed  a  diffuse  zone  with  a  moderate  bias  toward  the  metal. 
This  latter  result  provides  some  reassurance  that  small  simulations  can  give  usefully  information 
about  water  structure  within  two  or  three  layers  from  the  surface. 

Near  the  metal  the  ion  center  did  not  pass  z  =  0.4  nm.  This  meant  it  was  never  closer  than  0.3 
nm  to  the  point  where  the  wall  potential  passed  through  zero  (vertical  dashed  line  at  0.6  nm  in 
Figure  2).  At  this  point  of  closest  approach  there  is  still  room  for  one  water  molecule  between 
the  ion  and  the  metal.  In  Grahame's  picture  of  aqueous  electric  double  layers^’  ^  water  in  the 
primary  solvation  shell  of  small  ions  can  touch  the  electrode  surface.  In  contrast  in  the  model 
proposed  by  Bockris^’  ^  the  solvated  ion  does  not  displace  water  in  the  first  layer  next  to  tlie 
charged  metal,  and  the  ion  is  always  at  a  distance  of  about  two  water  molecules  from  the  surface. 
In  the  simulations  described  here  the  ion  behaved  like  the  picture  due  to  Grahame.  A  more 
detailed  model  of  adsorption  in  which  the  ion  interacts  with  atoms  and  electrons  comprising  the 
surface  is  beyound  the  scope  of  the  present  work.  In  passing  we  mention  that  metal-water  po- 
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tentials  exist  for  platinum  and  several  other  metals.  Generally  the  effective  well  depths  are  large 
and  to  calculate  distribution  profiles  requires  either  very  long  simulation  times  (>  10  ns)  and 
or  the  use  of  specialized  sampling  techiques^*^. 

The  water  profile  shows  some  new  structure  not  seen  in  water  without  ions^^’  Most  inter¬ 
esting  is  the  peak  closest  to  the  metal  surface  at  ca.  0.68  nm  due  to  a  few  localized  water  mol¬ 
ecules.  The  orientation  of  these  localized  water  molecules  can  be  determined  from  the  H_ST2 
and  PC_ST2  probability  distributions.  The  protons  on  these  molecules  give  rise  the  distinct  peak 
at  ca.  0.75  nm  in  the  H_ST2  distribution.  In  the  PC_ST2  distribution  there  is  peak  at  ca.  0.625 
nm  that  is  enhanced  more  than  the  second  peak  in  the  H_ST2  profile  at  ca.  0.6  run.  The  first 
water  peak  at  ca.  0.68  run  lies  between  the  first  H  and  PC  peaks  measured  firom  the  metal  sur¬ 
face.  The  positions  and  relative  intensities  of  these  peaks  suggests  that  some  of  the  localized 
water  molecules  have  one  proton  pointing  at  the  electrode.  This  would  permit  the  system  to 
adopt  an  ice-like  configuration  with  hydrogen  bonding  to  a  second  layer  of  water,  similar  to  that 
discussed  by  Lee  et  al*.  We  do  not  preclude  other  waters  with  two  protons  pointing  at  the 
surface  these  would  have  a  dipole  normal  to  the  surface  which  is  favored  on  the  basis  of  simple 
electrostatics. 

Figure  3  shows  the  electric  potential  (atom  xlO)  across  the  cell  calculated  by  the  atomic  charge 
method  with  a  test  function  with  gaussian  width  g  =  0.03  run.  Also  shown  are  the  components 
of  the  electric  potential  coming  from  the  ions  (monopoles)  and  the  point  dipoles  of  the  waters 
molecules.  The  component  potentials  were  calculated  by  the  molecular  method.  The  compo¬ 
nents  are  monopole  (mono),  dipole,  and  monopole  +  dipole  combined  (m+d  xlO).  Note  that  for 
Izl  <  0.6  run  the  monopole  and  dipole  terms  almost  cancel. 

The  monopole  contribution  (mono  in  Figure  3)  due  to  lithium  ion  alone  drops  monotonically 
as  expected  for  the  potential  inside  a  capacitor  where  the  charge  on  the  left  is  in  a  diffuse  layer 
spatially  separate  from  the  charge  on  the  right  plate  (metal)  at  z  =  0.931  nm.  Adding  the  dipole 
completely  changes  the  potential  (see  m+d  xlO  curve).  The  water  molecules  very  effectively 
screen  the  field  inside  the  capacitor,  except  for  the  region  z  >  0.68  nm  where  the  water  distrib¬ 
ution  drops  rapidly  to  zero.  The  peak  in  m+d  near  0.7  run  occurs  where  they  no  longer  cancel 
and  the  monopole  charge  dominates. 

Figure  4  shows  some  of  the  higher  order  electrostatic  components  of  the  electric  potential. 
Quadrupole  and  octopole  potentials  are  from  the  water  molecules.  To  allow  an  easy  comparison 
the  monopole  +  dipole  combined  (m+d)  potential  is  also  plotted.  Note  that  the  quadmpole  po¬ 
tential  is  approximately  constant  across  the  cell  except  near  the  metal  where  the  potential  has  a 
sharp  peak.  The  octopole  potential  is  everywhere  approximately  zero  except  near  the  metal 
surface.  Adding  the  quadrupole  (also  from  water  only)  term  brings  the  atomic  and  molecular 
calculations  into  some  measure  of  agreement  The  quadrupole  contribution  to  the  potential  is  a 
negative  constant  in  the  region  away  from  the  wall  potential.  This  occurs  because  its  contrib¬ 
ution  to  the  charge  <p„„;(r,/)>  contains  a  double  derivative  in  space  coordinates.  The  differences 
are  greatest  at  the  surfaces  where  the  atomic  method  traces  charge  density  smoothly  whereas  the 
molecular  method,  based  on  an  expansion  about  molecular  centers,  requires  many  high  multi- 
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poles  to  describe  the  field.  The  dangers  of  omitting  higher  raultipole  were  clearly  pointed  out 
by  Wilson,  Pororille  and  Pratt^^’  for  water  without  ions. 

We  close  this  section  with  a  discussion  of  changing  the  length  scale  over  which  the  test  function 
averaging  is  performed.  The  lithium  system  was  chosen  because  prior  worki*’  and  systematic 
electrochemical  experiments^  has  shown  that  small  ions  like  lithium  or  fluoride  form  diffuse 
screening  zones  near  charged  electrodes,  and  have  little  tendency  to  contact  adsorb.  In  Figure 
5  shows  the  relative  insensitivity  of  the  monopole  and  dipole  potentials  to  a  change  in  the 
gaussian  bin  widths  for  the  range  g  =  0.03  to  0.3  nm.  The  solid  curves  are  for  g  =  0.03  nm  and 
the  broken  curves  for  g  =  0.3  nm.  The  smearing  of  monopole  (ion)  charge  for  g  =  0.03  to  0.3 
nm  did  not  cause  it  to  overlap  the  metal  so  that  this  component  of  the  electric  potential  hardly 
changed.  The  dipole  potential  changed  most  with  g,  notably  for  z  >  0.68  nm  since  the  water 
density  extended  all  the  way  to  the  repulsive  region  of  the  wall  potential  of  the  metal.  Since  the 
monopole  and  dipole  potential  almost  cancel  in  the  solution  phase  this  makes  the  m+d  potential 
sensitive  to  the  detailed  behaviour  of  the  dipole  near  the  surface. 

Figure  6  shows  the  sensitivity  of  higher  order  electrostatic  multipole  potentials  to  gaussian  bin 
widths  for  the  range  g  =  0.03  to  0.3  nm.  Solid  curve  for  g  =  0.03  nm,  long  broken  curve  for  g 
=  0.1  nm,  and  short  broken  curve  for  g  =  0.3  nm.  The  surface  features  of  the  quadrupole  and 
octopole  are  washed  out  by  the  value  g  =  0.3  nm  the  size  of  the  water  molecule.  Note  that  the 
quadrupole  potential  is  averaged  to  a  finite  value,  very  roughly  constant  over  the  cell.  The  po¬ 
tential  from  the  octopoles  averages  to  zero  over  the  whole  cell  for  bin  widths  greater  than  g  — 
0.1  nm. 

The  effect  of  changing  the  gaussian  width  on  the  total  potential  is  shown  in  Figures  7  and  8. 
In  Figure  7  the  sensitivity  of  the  atom  potential  to  the  gaussian  bin  width  is  shown  (solid  curves) 
for  g  =  0.03  to  0.3  nm.  The  broken  curves  show  the  corresponding  monopole  potentials.  The 
difference  between  the  two  sets  of  curves  is  due  to  the  water.  Note  that  the  peak  in  the  potential 
near  0.8  nm  is  completely  washed  out  at  g  =  0.3  nm.  Figure  8  shows  a  detail  of  Figure  7  in  the 
range  of  small  electric  potentials  (-.05  to  .05).  The  gaussian  bin  widths  are  the  same  as  before 
g  =  0.03  to  0.3  nm.  The  corresponding  monopole  potentials  are  plotted  for  reference.  Note  that 
the  peak  in  the  potential  near  0.75  nm  is  completely  washed  out  by  g  =  0.3  nm. 

Figures  7  and  8  demonstrate  that  as  the  size  of  the  bin  is  increased  from  g  =  0.03  nm  (.1  x  water 
molecule  dimension)  to  g  =  1.0  nm  (  3  x  water  molecule  dimension)  the  structure  in  the  electric 
potential  near  the  surface  due  to  the  water  profile  is  lost.  The  total  potential  becomes  monotonic 
and  resembles  the  shape  of  tlie  monopole  ptotential  curves  calculated  using  the  monopole  charge 
distribution.  Since  the  component  of  the  water  dipole  perpendicular  to  the  surface  is  not  aver¬ 
aged  to  zero  tlie  atom  potential  at  the  surface  always  remains  smaller  in  magnitude  than  me 
monopole  potential.  We  can  regard  this  as  an  expression  of  dipole  dielectric  polarization  which 
remains  after  gaussian  averaging  and  reduces  the  value  of  the  surface  potential  by  roughly  an 
order  of  magnitude  (curves  for  g  -  .3  nm).  As  expected  the  monopole  potentials  are  not  as 
sensitive  to  the  value  of  the  width  g  as  the  atomic  potentials.  This  is  an  interesting  qualitative 
result  because  it  shows  the  transition  from  the  microscopic  scale,  where  surface  water  oscillatory 
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structure  dominates  the  potential,  to  macroscopic  scale  behaviour  where  water  contributes  simple 
scaling  of  the  electric  potential. 
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IV.  CONCLUSIONS 

For  a  system  with  non  contact  adsorbing  cation  we  have  used  two  methods  to  calculate  a  time 
independent  electric  potential  from  atom  trajectories  generated  by  a  constant  (N,V,T)  molecular 
dynamics  simulation.  A  key  feature  of  the  simulation  was  the  calculation  of  all  long  range 
coulomb  interactions  including  all  electric  image  interactions  without  truncation.  Each  method 
provides  some  insight  into  the  components  of  the  potential  contributed  by  the  ion  and  the  water 
molecules  and  the  sensitivity  of  the  multipole  fields  to  the  scale  over  which  spatial  averages 
were  performed.  The  variation  in  potential  near  the  metal  was  due  to  partially  oriented  water 
molecules.  The  importance  of  water  dipole,  quadmpole  and  octopole  potentials  was  highlighted. 
Quadrupoles  contribute  roughly  a  constant  value  to  the  potential  in  the  bulk  and  oscillatory 
stmcture  near  the  surface.  A  significant  contribution  from  octopoles  and  higher  order  multipoles 
occured  near  the  surface  where  the  fields  are  strong  and  the  potential  rapidly  varying.  The  in¬ 
sensitivity  of  monopole  and  dipole  to  the  spatial  bin  width  was  in  contrast  to  the  behaviour  of 
the  quadmpole  and  octopole. 
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Figure  1.  Schematic  diagram  of  the  charged  metal  electrode  aqueous  electrolyte  interface.  The 
metal  is  shown  as  a  flat  negatively  charged  surface  and  there  is  a  fully  hydrated  cation  at  the 
outer  Helmholtz  plane  (OHP)  two  water  molecules  distant.  The  top  part  shows  a  schematic  of 
the  electric  potential  assuming  the  ions  are  point  charges  in  a  dielectric  continuum  restricted  to 
the  solution  side  of  the  OHP. 
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Figure  3.  Electric  potential  (atom  xlO)  and  components  of  the  potential  calculated  by  the  mo¬ 
lecular  method.  The  components  are  monopole  (mono),  dipole,  and  monopole  +  dipole  combined 
(m-i-d  xlO).  Note  that  for  Izl  <  0.6  nm  the  monopolc  and  dipole  tenns  almost  cancel.  Simulation 
cell  contains  one  Li-"  ion  and  157  st2  waters  with  the  metal  electrode  on  right  hand  side,  and 
dielectric  on  the  left.  Image  plane  at  z  =  0.931  nm  and  wall  potentials  go  through  zero  at 
lzl=0.682  nm. 
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Figure  4.  Higher  order  electrostatic  components  of  tlie  electric  potential.  Quadrupole  and 
octopole  potentials  from  the  water  molecules.  For  comparison  die  monopole  +  dipole  combined 
(m+d)  potential  is  also  plotted.  Note  that  the  quadrupole  potential  is  approximately  constant 
across  the  cell  except  near  tlie  metal  where  the  potential  has  a  sharp  peak.  The  octopole  potential 
is  everywhere  approximately  zero  except  near  the  metal  surface.  Simulation  cell  contains  one 
Li"^  ion  and  157  st2  waters  with  the  metal  electrode  on  right  hand  side,  and  dielectric  on  the  left. 
Image  plane  at  z  =  0.931  mn  and  wall  potentials  go  through  zero  at  lzi=0.682  mn. 
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Figure  5.  Insensitivity  of  low  order  electrostatic  multipole  potentials  to  gaussian  bin  widths  for 
the  range  g  =  0.03  to  0.3  nm.  Solid  line  result  for  g  =  0.03  nm.  broken  line  result  for  g  =  0.3 
nm.  Monopole  and  dipole  do  not  change  value  except  the  dipole  for  z  >  0.68  nm  where  the 
water  density  is  very  small.  Simulation  cell  contains  one  Li"'^  ion  and  157  st2  waters  with  the 
metal  electrode  on  right  hand  side,  and  dielectric  on  the  left.  Image  plane  at  z  =  0.931  nm  and 
wall  potentials  go  through  zero  at  lzl=0.682  nm. 
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Figure  6.  Sensitivity  of  higher  order  electrostatic  multipole  potentials  to  gaussian  bin  widths  for 
the  range  g  =  0.03  to  0.3  nm.  Note  that  the  quadrupole  loses  the  surface  peak  for  a  width  equal 
to  the  water  molecule  dimension  0.3  nm.  The  octopole  potential  averages  to  zero  over  the  whole 
cell.  Simulation  cell  contains  one  Li^  ion  and  157  st2  waters  with  the  metal  electrode  on  right 
hand  side,  and  dielectric  on  the  left.  Image  plane  at  z  :=  0.931  nm  and  wall  potentials  go  through 
zero  at  lzl=0.682  mn. 
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Figure  7.  Sensitivity  of  the  atom  method  electric  potential  to  to  gaussian  bin  widths  for  the  range 
g  =  0.03  to  0.3  nm.  The  corresponding  monopole  potentials  are  plotted  for  reference.  Note  that 
the  peak  in  the  potential  near  0.8  nm  is  completely  washed  out  at  g  =  0.3  nm.  Simulation  cell 
contains  one  Li+  ion  and  157  st2  waters  with  the  metal  electrode  on  right  hand  side,  and 
dielectric  on  the  left.  Image  plane  at  z  =  0.931  nm  and  wall  potentials  go  through  zero  at 
izl==0.682  nm. 
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Figure  8.  Detail  showing  the  small  potential  range  (<  .05)  of  tiie  sensitivity  of  tlie  atom  metliod 
electric  potential  to  to  gaussian  bin  widths  for  the  range  g  =  0.03  to  0.3  nm.  The  corresponding 
monopole  potentials  are  plotted  for  reference.  Note  tliat  the  peak  in  tlie  potential  near  0.8  nm 
is  completely  washed  out  at  g  =  0.3  nm.  Simulation  cell  contains  one  Li^  ion  and  157  st2  waters 
with  the  metal  electrode  on  right  hand  side,  and  dielectric  on  the  left.  Image  plane  at  z  =  0.931 
nm  and  wall  potentials  go  through  zero  at  lzl=0.682  nm. 
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